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Abstract 



It has been speculated that the recently detected ultra-high energy cosmic 
rays may originate from the decays of relic particles with mass of order 10^^ 
GeV clustered in the halo of our Galaxy. This hypothesis can be tested 
through forthcoming measurements of the spectra of both high energy cosmic 
nucleons and neutrinos, which are determined in this model by the physics 
of QCD fragmentation, with no astrophysical uncertainties. We evolve frag- 
mentation spectra measured at LEP energies up to the scale of the decaying 
particle mass by numerical solution of the DGLAP equations. This enables 
incorporation of the effects of supersymmetry on the development of the cas- 
cade and we also allow for decays into many-particle states. The calculated 
spectral shape agrees well with present cosmic ray data beyond the Greisen- 
Zatsepin-Kuzmin energy. 
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I. INTRODUCTION 



Within the last decade, the Fly's Eye ^ and its successor HiRes Q, as well as the Akeno Giant 
Air Shower Array (AG AS A) j^, have detected over 70 ultra-high energy cosmic rays (UHECRs) 
with energies measured reliably to be in excess of Eqzk — 4 x 10^^ eV — the Greisen-Zatsepin- 
Kuzmin (GZK) 'cutoff' energy |^. The older air shower arrays at Volcano Ranch |p, Haverah Park 
1^], and Yakutsk have added another 45 (recalibrated) post-GZK events. At energies above 
Eqzk, the typical range of protons falls rapidly because of interactions with cosmic microwave 
background photons ^]^, becoming as low as ~ 20 Mpc at 3.2 ± 0.9 x lO^o eV, the highest energy 
event observed by Fly's Eye; the range is even smaller for heavy nuclei |jlO|.Q For photons with 
energy above Eqzk , the dominant opacity is due to scattering on the extragalactic radio background 
which is poorly known experimentally |16|; the mean free path is estimated to be ~ 1 — 5 Mpc at 
10^^ eV [^J^ There is no restrictive constraint on the range if the UHECR primaries are neutral 
and weakly interacting, viz. neutrinos which have only a tiny probability for scattering on the 
relic cosmic neutrino background [ p^j2^ f| This also means however that neutrinos cannot initiate 
the observed airshowers as their cross-section for deep inelastic scattering on nucleons is far too 
small |27].0 Thus all the available evidence suggests that cosmic rays with energies E > Eqzk are 
protons, in which case they must originate within the Local Supercluster.0 



^The shower elongation rate measured by Fly's Eye [|i|,[TT|| and HiRes indicates a change in 
composition from heavier nuclei towards nucleons at the highest energies. Muon data from AG AS A 



1 13 1 are consistent with this trend when analysed using the same hadronic interaction model |14,15| 



^The highest energy event Q had a depth of maximum of 8I5I35 gem ^ cf. the expected value 



of 1075 gem for a photon primary |18|. The muon content of horizontal showers in the Haverah 



Park data implies an upper bound of 65% on the photon component, and a (hadronic interaction 



model dependent) upper bound of ~ 50% on the iron nuclei content, of post-GZK UHECRs |19|. 
It should be emphasised however that a complete understanding of air shower development is still 
lacking, so the composition of UHECRs is not definitively established yet [pC|j2 



^However if the background neutrinos have mass of ~ 0.1 — 1 eV, then UHE cosmic neutrinos 
can annihilate resonantly on these to generate 'Z-bursts' which create photons and nucleons above 
Eqzk [p^- To match the observed UHECR flux given experimental limits on the UHE cosmic 



neutrino flux |25| requires however a substantial increase in the local density of relic neutrinos |26|. 



^Neutrinos may interact more strongly if dramatic new physics is invoked e.g. extra dimensions 
and quantum gravity at the TeV scale [28|. However explicit calculations indicate so far that 
the increased cross-section would still be inadequate to create the observed air-showers |p9|. This 



possibility is testable through studies of penetrating airshowers [30 1 and upward-going muons |31|. 



^A suggested correlation between UHECR arrival directions and cosmologically distant radio 
sources [^] motivates the possibility that the primary is a new heavy hadron (e.g. a super symmetric 



uds-gluino state) with a higher GZK cutoff |33|. The correlation has however been questioned |34|; 



moreover there are stringent experimental bounds on new stable strongly interacting particles pi 
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However they cannot originate in the Galactic disk since significant anisotropies would then be 
expected at post-GZK energies |^^, while the observed sky distribution is consistent with isotropy 
A number of doublet and triplet events, collimated within the experimental angular 
resolution of a few degrees, have been observed but the chance probability of these arising from an 
isotropic distribution is estimated to be of order 10% ||4^, hence statistically not significant|] 

There have been a plethora of suggestions as to the possible origin of UHECRs, which may be 
broadly classified into 'top-down' or 'bottom-up' models. In the former class, the extreme energies 
are supposedly provided by the decay of relic topological defects (TD) |4l[-43,15 or super-massive 
particles [46-49|. In the latter class, the UHECRs are assumed to have been accelerated up from low 
energies in astrophysical sites such as gamma-ray bursts [50|, active galaxies [^] or super- massive 
black holes in dormant quasars 15^. In fact there are severe difficulties with accelerating particles 



to the highest observed energies by any known astrophysical mechanism |53|. Moreover even if such 
sources are distributed nearby like the observed galaxies, the injection spectrum must be harder 
than dA^ /dE cx E~'^ in order to overcome the GZK energy losses and match the measured UHECR 
flux [p|,|5^, so the energetic demands are considerable. Gamma-ray bursters may possibly have 
the required energy but, being more abundant at high redshift, cannot yield the locally observed 
UHECR flux |^5|. There are a few active galaxies within the 'GZK sphere' but even if they can 
accelerate UHECR these would need to be deflected through rather large angles (and isotropised) 
by the intergalactic magnetic field [^6|; this requires the field strength to be ~ 10^ times stronger 
than the upper limit of order nanogauss usually inferred from observational bounds on Faraday 
rotation in distant radio sources |^7|. There is ongoing discussion concerning the possibility of such 
high fields 



as well as of specific sources such as Cen A |59|-|61| and Virgo A [|62|,|63| . Although 
there is no compelling astrophysical model for post-GZK UHECRs, some of these suggestions are 
interesting in their own right and can be tested with forthcoming data. 

In this paper, we investigate the spectral signature of 'top-down' models. Relic topological 
defects will in general be cosmologically distant and release particles at the GUT-scale, resulting in 
a typically excessive fiux of secondary low energy 7-rays as the primaries interact with intergalactic 
radiation fields ^,45|. So we favour the possibility that UHECRs originate from the decay of super- 
massive relic particles which, being cold dark matter (CDM), are locally clustered in the halo of 
our Galaxy, thus naturally evading the GZK-cutoff and also presenting an approximately isotropic 
distribution |48,49|. Such particles must have a mass ^ 10^^ GeV to account for the highest energy 
event, and a lifetime exceeding the age of the universe. A physically well-motivated candidate with 
both the required mass and lifetime is the "crypton" — the analogue of a hadron in the hidden 
sector of supersymmetry breaking in string theories [341.^ It has been noted that such particles 
can readily be produced with a cosmologically interesting abundance through the time-varying 
gravitational field at the end of infiation [66|, or possibly, during (re)heating after inflation |67|. 



^At around 10^^ eV Fly's Eye ^ and AGASA have detected anisotropy towards the Galactic 
plane, consistent with the usual belief that cosmic rays at these energies have a galactic origin. 

^However if only the paired events within ±10" of the supergalactic plane are considered, the 
chance probability that they arise from an isotropic distribution is less than 1% |4C 



^Other particle candidates which may be sufficiently long-lived have also been proposed [S5|. 
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A key test of this model is the expected smah anisotropy due to the offset of the Sun from the 
centre of the Galaxy |68|.0 Some authors have argued |^,^ that present limits on the anisotropy 



already place severe constraints on the model; however a recent detailed analysis of the expected 
anisotropy signal ||7l| of post-GZK UHECRs shows it to be both compatible with all present data 
and moreover capable of being definitively detected by the Pierre Auger Project which is presently 
under construction ||7^, and by the proposed Extreme Universe Space Observatory These 
experiments will provide a very substantial increase in the presently meagre statistics of post-GZK 
UHECRs and also be able to distinguish between incident photon, neutrino and nucleon primaries. 
This is clearly an opportune time to make predictions of the expected spectrum and composition 
of UHECRs in the decaying dark matter (DDM) model. 

II. THE DECAYING DARK MATTER MODEL 

The essential idea ||^,49| here is that because of gravitational clustering, very massive relic 
particles would naturally have a density in the halo of our Galaxy which is enhanced by a factor 
of ~ 10^ over the cosmic average. This is approximately the same factor by which the distance to 
the horizon exceeds the halo radius, so if these particles decay to produce UHECRs the flux from 
cosmologically distant regions will be at best comparable to the flux from our halo even if there 



is no absorption |74]. Given the GZK attenuation, the extragalactic flux of protons and photons 
in this model would be quite negligible compared to the flux from the Galactic halo, whereas for 
neutrinos there would be an extragalactic component comparable to the galactic one [23| ]. 

The injection spectrum produced by the decay of a population of superheavy dark matter 
particles X, with number density nx, is proportional to the inclusive decay width of X: 

Tx Tx dE 

where we have taken the lifetime of X to be longer than the age of the universe, tx ^ ~ 10^'' yr. 
Photons produced in the galactic halo can interact with the galactic low frequency radio background 
on their way to the Earth; the number of high energy photons reaching the Earth may then be 



significantly reduced. However the radio background is presently very poorly known [16| so we do 
not include photon processing in the halo in this work. We consider a spherical halo of radius i?haio 
and uniform density njf'° (see Ref. |71] for a thorough study of UHECR and halo models). The 



galactic halo contribution to the UHECR spectrum is then given by 

_^halo(^) = ^i?halo^""'°(i?). (2) 

In order to calculate the inclusive decay width, some assumptions have to be made about the 
microphysics of X decay. It is usually assumed in top-down models (i.e. for both TD and DDM) 
that the initial decay is into a parton-antiparton pair. Any angular dependence of the decay 
amplitude is immaterial since the flux at Earth is close to isotropic. Therefore in the following 



Section (IH) we shall consider decays into quarks and gluons (and also their superpartners when 



allowing for the presence of supersymmetry). 



^If the clustering in the arrival directions turns out to be statistically significant, it may be an 



indication of dumpiness in the halo dark matter distribution, as is indeed expected for CDM |69| 
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III. EVOLUTION OF FRAGMENTATION FUNCTIONS 



A. Inclusive Decay Width 

For a particle X with mass Mx , decaying into partons whicli hadronise into particles of type h 
(carrying a fraction x of the maximum available momentum Mx/2, and a fraction z of the parton 
momentum), the inclusive decay width can be factorised as 

1 dr{X^h + ...) f'dzj_ dra{y,^l^,M^] 

dx „Jx z Ta dy 



D^ai^.l^'). (3) 

y=x/ z 



The first factor, the decay width of X into parton a, dFa/dy, is calculable in perturbation theory; 
in lowest order and for 2-body decay it is proportional to 5{1 — y). The second factor, the non- 
perturbative D^, is the fragmentation function (FF) for particles of type h from partons of type a. 
Both factors depend on the factorisation scale The dependence on /i of the physical inclusive 
decay width would cancel out if the calculation could be carried out to all orders in perturbation 
theory but in a finite-order calculation some dependence on fi remains to orders higher than those 
calculated. The choice /i ~ Mx is the most appropriate |77|. Every perturbative term in the above 



factorisation depends also on an arbitrary renormalization scale; following the usual convention, 
we take the renormalization scale to equal the factorisation scale. 

The cancellation of the dependence on /i between the perturbative factors and the FFs con- 
strains the latter to satisfy the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations 



1 78, 79 1 . Given experimental data on hadronisation at some low energy scale, say Mz, an initial 
set of fragmentation functions D^{x, M"^) can be extracted. These can then be evolved using the 
DGLAP equations to obtain the FFs at some higher scale D^{x, M'j^). This evolution gives rise to 
scaling violations in the x-dependence of the inclusive decay width of X in exactly the same way 
that scale-dependence of structure functions and FFs produce scaling violations in experimentally 



measured hadronic cross-sections |76 

The DGLAP equations can be written as 

5^1^ = E ^na(x,a.(/x2)) ® Di{x,f,% (4) 
n /X ^ vr 

where Q;s(/i^) is the strong coupling constant and Pba{x, Os) is the splitting function for the parton 
branching a ^ b. Here the convolution of two functions A{x) and B{x) is defined as 

A{x)®Bix)= / ^A{z)B{^). (5) 
The splitting functions can be expanded perturbatively: 

Pbaix,as) = Pbaix)+0{as). (6) 

We limit our study to leading order in Os and therefore ignore 0{as) corrections to the splitting 
functions. It is also convenient to define the following dimensionless evolution parameter 

b being the coefficient in the leading order /3-function governing the running of the strong coupling: 
/3{as) = —ba^. We take to represent the sum of particle h and, if different, its antiparticle h. 
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B. Standard Model Equations 



The Standard Model DGLAP equations for the evolution of fragmentation functions are well- 
known 1^6],^]. There are two parton species: quarks qk, k = I, . . . np, and gluons g, with np the 
total number of flavours. Conventionally one defines the following linear combinations (for ease of 
notation the superscript h is omitted): 



D„ = 



k 

D. 



Qk 



The non-singlet functions D - and Dq^, obey the equations 



(8) 
(9) 
(10) 

(11) 



Pr 



qq 



Pr 



qq 



(12) 
(13) 



while the evolution of the singlet function Dq is coupled to that of the gluon function Dg as 



dr 



Pqq 2nFP< 



gq 



P 



qq 



P 



qq 



Dq 



(14) 



The splitting functions were calculated in Refs. [78,7^. Given the FFs at some initial scale /^o for 
the quarks and gluon g, we can now determine their evolved values at some other scale fi to 
leading order in as, using Eqs. (12-14). 



C. Supersymmetric Equations 

In a supersymmetric (SUSY) model, besides the quarks and gluons one has their superpartners: 
squarks Sk and gluinos A. In addition to the linear combinations (^-pT]) we now define 

D.^Ds, + D,^, (15) 

k 

Ds^T^Dsp (16) 

k 

D^-^Ds,-D,^, (17) 
Ds, = - —Ds. (18) 
The non-singlet function D^- and D^- evolve together, as do Dg^. and -Ds^.: 
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The singlet functions for quarks and squarks, Dq and Ds, are coupled to the gluon and gluino 



functions, Dg and Dx, as 



fDq\ 

Ds 
\DxJ 



Pr 



qq 



qg 

\PqX 



2nFPgq 

p 

2npPgs 



P, 



-I sq 
P 

J^sg 
P 

^ ss 
PsX 



2nFPxq \ 

Pxg 
In^Pxs 

Pxx J 



(Dq\ 

Dg 
Ds 
\DxJ 



(21) 



In leading order, Eqs. (19-21) allow us to calculate the fragmentation functions for all quark and 
squark flavours, gluons and gluinos at some scale /i, given their values at some initial scale /io- 



Incidentally, using the following relationships among the SUSY splitting functions [S0|, 



Pqq ~l~ Psq 

Pgq + Pxq 
Pqg + Psg 

^99 + -^^3 



Pss ~l~ Pqs^ 
Pxs ~l" Pgs, 
PsX + PqX, 
Pxx + PgX, 



one can reduce the 4x4 matrix equation (|2^) to the 2x2 matrix equation 



9. 



Dq 
Dn 



Pqq Pqs 2np(^Pgq P^ 



gs) 



Pr 



99 



P 



qX 



P. 



99 



P 



9A 



Dq 
Dn 



(22) 
(23) 
(24) 
(25) 



(26) 



The SUSY DGLAP equations have been given in the literature to leading order for structure 
functions pO| , p| ]. Here we have presented their form for FFs. It is easy to see that (except for 
Eq. (|26|)) one just needs to transpose the matrix elements, keeping the n-p factors in the same place, 
to move from structure function to fragmentation function equations. For Eq. (p6[), the relative 
minus sign between q, s and g, A changes to a plus sign and there is a rearrangement of the splitting 
functions in the 2x2 matrix. The SUSY splitting functions were calculated in Refs. 80|,8l[. 



D. Numerical Evolution: Algorithm and Initial Conditions 

Several numerical algorithms and their implementation codes are available to solve the set of 
integro-differential DGLAP equations for structure functions when quarks and gluons are the only 
particles inside the hadrons ||8^]8^ . The super symmetric evolution of structure functions has been 
considered in Ref. |84|. These codes are suitable for collider physics but are rather unwieldy to 



use for studying cosmic ray production. We have generalised the Laguerre method |82| to include 
SUSY evolution for fragmentation functions and written a numerical code to calculate cosmic ray 
production from superheavy particle decay. The main idea of the Laguerre method is to expand 
Eq. (^) in Laguerre polynomials; the coefficients of the series expansion are given by simple recursive 
relations appropriate for computer calculations. These recursive relations have been given so far 
for the 1-dimensional case (Eqs. ^-13) and the 2-dimensional case (Eq. [l^ ). However to study the 
general SUSY evolution requires one to work out the recursive relations for the 4-dimensional case 
(Eq. |2^). This is not a trivial step if one follows the procedure given for lower dimensional cases 
in the literature. Hence we have calculated the general recursive relations for any dimension. The 



details of this procedure, as well as our numerical code, will be given elsewhere |85|. 

We begin with FFs at the energy scale /xq = Mz and evolve them to the final energy scale 
/i = Mx using Eq. (^). For the initial fragmentation function of baryons, DP{x, M"^) + D2{x, M"^), 



7 



we adopt the fit performed in Ref. |86| to LEP hadronic data as sliown in Fig. |l[ For photons 
and neutrinos we generate initial data at the Z peak using the QCD Monte Carlo event generator 
HERWIG |89|. Comparison with LEP data shows that although HERWIG overproduces baryons 



at high X [|86| , p7[ , its photon and meson output at the Z peak matches the experimental spectra 
remarkably well. Since neutrinos mainly come from charged pion and kaon decays, one can thus 
be confident in taking the HERWIG generated FF for neutrinos as the initial condition for the 
evolution. There is also a sizable contribution from heavy flavour decays to the neutrino spectrum 
at high X which is explicitly taken into account by HERWIG. Figure |2| shows our initial spectra at 
the energy scale Mz for photons and neutrinos together with LEP data for photons. 

We assume flavour universality in the decay of X, hence we consider only the coupled singlet 
quark and gluon evolution equations (|lj) for the Standard Model (SM), and the coupled singlet 
quark, gluon, singlet squark and gluino evolution equations ( ^l] ) for a SUSY model. A (s)parton is 
not included in the evolution as long as the energy scale is lower than its mass; when the threshold 
for its production is crossed, it is added to the evolution equations with an initially vanishing FF 
and it is assumed to be a relativistic particle. 

In the SM case we evolve the q and g initial fragmentation functions from Mz to Mj, the top 
quark mass, with the number of flavours set to np = 5, and then evolve from Mt to Mx with 
np = 6. (However assuming np = 6 in the whole range from Mz to Mx does not introduce any 
significant difference in the final spectrum.) 

In the SUSY case we evolve the q and g initial fragmentation functions from Mz to the su- 
persymmetry breaking scale MsusY > Mt using the SM equations to obtain D'I {x , Mgugy) , with 
i = q,g. Then we take D^{x, M^^^y), i = q,g, and Dj{x, Mg^gy) = 0) J = s,X, and evolve them 
from MsusY to Mx using the SUSY equations. All spartons are taken to be degenerate with a 
common mass MgusY- (In the context of structure functions, a SUSY model with different masses 
for s and A was studied finding no significant difference with models having the same mass 
for all superpartners. One might expect the same result to hold for FFs.) 

For very small x values, x ^ 0.001, coherent gluon emission becomes important. In this 
regime Eq. (^) no longer holds and the kernel of the integral must be modified as Da{x/z,ii^) — >■ 
Da{x/ z'^fj,'^). This is the Modified Leading Logarithm Approximation (MLLA) |7^,9^ which has 
been used by other authors to calculate the cosmic ray spectrum in top-down models p3,48,91 



However the MLLA spectrum cannot be normalised since a very small (and uncertain) fraction of 
the total energy is released in this kinematic region; moreover in the context of the DDM model for 
a particle mass of order the hidden-sector scale, the observed highest energy cosmic rays correspond 
to large values of x |4£,92]. There are also other corrections to the DGLAP equations arising from 
a more rigorous treatment of threshold effects |77] and, of course, next-to-leading order corrections 
to as which we have not considered. Nevertheless, the level of precision afforded by the DGLAP 
equations in leading order, when calculating scaling violations in the inclusive decay of X, is quite 
adequate for comparison with the present experimental data on UHECRs. 

For X ~ 1 there are large uncertainties (see Fig. |^) in the experimental determination of the FFs 
|l88|| . Generating values for these using HERWIG requires rather large amounts of computational 
time since very few particles fall in the bins with x > 0.8, if the number of simulated events is kept 
to a reasonable number |^9|. Furthermore, the algorithms which solve the DGLAP equations show 
poor convergence in this regime. These problems motivate an alternative approach to calculating 
high X fragmentation. One can take a power-law fit to the experimental data at AIz, and then use 
the Mellin transform technique |76| to calculate its evolved value at Mx- If at tq one has the fit 
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Da{x, To) = K{1 — x)^, where K and k are constants, then when x — > 1 we find 



Da{x,T)=K 

where in the SM 

A, 



Qs(to) 



Ba 



T{k + l) 



T{Aa{T-To) + k + l) 



(1-x) 



A, 



and in a SUSY model 



Aq — Ag 
Aa = A^ 



2Cf, 
2Ca, 

= 2Cf, 
= 2Ca, 



R — ct 

~ 2-Kh 



27 



1 

27rb 



Ca - 27) - fTRnp 



Bq = B\ 



1 

27rfe 



Ca (I - 27) - Trtif 



(27) 

(28) 
(29) 



(30) 
(31) 



The 5C/(3) fundamental and adjoint group factors are respectively Cp = 4/3, Ca = 3, 6 is the 
leading order coefficient of the /5-function (|^, Tr = 1/2, and 7 = 0.57721 ... is the Euler constant. 
Eq. ( |27| ) shows that for x ^ 1, a power-law evolves at higher energy into a power-law with a 
bigger power index: the higher the energy the faster the drop in (1 — x). 



E. Numerical Evolution: Results and Discussion 

First we present the evolution of the baryon fragmentation function in the SM. In Fig. ^ we 
show the fragmentation functions for baryons {p and n) from quarks and gluons at the scale Mz 
as measured at LEP, and their evolved shape at Mx = 10^*^, lO^'^ and 10^^ GeV. (Following the 
standard convention we always plot the quantity x^Da{x, /j,'^), which is proportional to the cosmic 
ray energy spectrum E^J{E).) As the final scale increases, the number of particles grows at low x 
and decreases at high x as is well-known from many previous studies of scaling violations. 

Next we compare the SUSY evolution of FFs with SM evolution. In Fig. ^ we show the 
common initial baryon FF at Mz, its shape after SM evolution up to Mx = 10^^ GeV, and its 
evolved shape at the same final scale after SUSY has been switched on at MsusY = 400 GeV. It 
is clear that the SUSY curve has evolved further than the SM curve, chiefly due to the different 
running of as{fi'^). In a SUSY model Og decreases with increasing energy more slowly than in 
the SM because of the increased contribution to the /5-function from the SUSY partners. Since 
the rate of change di^^2Da{x, fj?) is proportional to (see Eq. Q), a bigger as translates into 
a larger amount of evolution. In other words, given the same initial and flnal scales we obtain 
tsvsy{Mx) > tsm{Mx), using Eq. (0). 

Figure ^ shows the quark and gluon FFs at Mz, their evolved values at Mx = 10^^ GeV using 
the SUSY equations for scales larger than MgusY , and the radiatively generated squark and gluino 
functions, all at the same final scale Mx- We find that starting from vanishing values at MgugY 
the squark and gluino functions start to grow and catch up rapidly with the quark and gluon 
functions, respectively, at small x. This behaviour can be understood qualitatively if one bears 
in mind that at low x the leading splitting function for quarks is 2nYPgq ~ ^u-pCy/x, which is 
equal to the leading splitting function for squarks 2n-pPgs ~ AtiyC-p/x. For gluons and gluinos the 
leading splitting functions tend as well to a common value, Pgg ~ 2Ca/x and Pgx ~ 2Ca/x, which 
is however different from that of quarks and squarks. 
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SUSY evolution does not depend strongly on the chosen super symmetry breaking scale MsusY- 
In Fig. ^ we show the curves obtained taking MsusY = 200, 400 GeV and 1 TeV. The higher 
the value of MsusYi the less evolved the final curves for q and g. This follows from our earlier 
comparison of SM vs. SUSY evolution. If SUSY switches on later (higher Msusy) the energy 
range over which the SM equations hold is larger. As we have seen already, DGLAP evolution is 
slower when just the SM equations are employed. 

The evolution of the photon and neutrino spectra can also be studied using the DGLAP equa- 
tions. We take HERWIG simulated data for the initial FFs for photons and neutrinos (sum over 
all three flavours) at Mz and use the SM and SUSY evolution equations ( p^ ) and ( [2l| ) to obtain 
their spectra at Mx- We plot in Fig. ^ the fragmentation functions for baryons, photons and 
neutrinos for a decaying particle of mass Mx = 10^2 Qgv, in both the SM and in a SUSY model 
with Msusy = 400 GeV. As pointed out earlier we assume flavour symmetry in the decay of X; 
moreover all the calculated FFs are colour-weighted. Hence for the SM the total contribution to 
the FF for particle of type h is just the sum of quark singlet and gluon FFs 

D^ = D^ + D^g, (32) 

while in a SUSY model we take the sum of quark and squark singlets, gluon and gluino FFs 

= + D^ + + Di . (33) 



IV. ULTRA HIGH ENERGY COSMIC RAY SPECTRUM 

We can now translate the calculated fragmentation functions into the expected cosmic ray spec- 
trum in order to confront the observational data. In the previous Section (|IIID we have calculated 
quark singlet and gluon functions for the SM, and quark singlet, gluon, squark singlet and gluino 
functions for SUSY. In the absence of a specific model for the different branching ratios we have 
weighted all the (s)parton contributions evenly but for colour weights as in Eqs. (^^33). Thus the 



quark contribution is the most important as seen from Figs. ^ and ^; in the SUSY case the squark 
contribution is also relevant for x < 0.01. For 2-body decay, x = 2E/Mx-, so from Eqs. (|l|-|3|) we 
have the following expression for the flux of particle h: 

E3j^^^°{E) = Bx^D^{x,Mj^). (34) 

We have multiplied the flux by -E^, as is usual, to emphasise the structure in the spectrum near the 
GZK energy. The normalisation factor B is common for the galactic halo flux of baryons, neutrinos 
and photons, and determines the quantity nx/rx (see Eq. ^). 

Let us now compare the calculated cosmic ray flux to the published data from Fly's Eye 0, 
AGASA Haverah Park and Yakutsk Q. In Ref. |75] these data have been carefully assessed 



for mutual consistency and appropriate adjustments made to the energy calibration, e.g. the 
original AGASA data [0] was reduced in energy by 10% to match the Akeno 1-km^ array data 



which covers the better explored energy region ~ 10^^ — 10^^ eV [93|. The authors recommend 
adoption of the following standard differential energy spectrum below the GZK energy, in the 
range 4 x 10^"^ eV < S < 6.3 x 10^^ eV: 

/ E \ -3.20±0.05 

JiE) = (9.23 ± 0.65) X lO-^^ m-^s-^sr-ieV-^ f _____ j , (35) 
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with the spectrum flattening at higher energies as J{E) oc E 2-75±o.2 ^p^Q |^]^g GZK energy, 
and extending further to at least 3 x lO^'^eV |73]. Thus the UHECR spectrum can naturally 



be interpreted ^ as the superposition of the 'low energy' component (35), and the new 'flat' 
component that extends into the post-GZK region. The former is presumably galactic in origin 
(consistent with the detection of anisotropy at ~ 10^^ eV [ p7| , pS| ]), while the latter is interpreted 
|p9|| as produced by the decay of a superheavy particle population in the galactic halo. Taking 
baryons to be the dominant primary UHECRs as indicated by experiment, the total flux is 

E'^J{E) = ^ + Bx^D^^'y^'^ix, Ml), (36) 

where the values of k and m can be read off Eq. (^). Note that since /^^aryon j^^y j^g^yg 
similar shape, taking photons to be the primaries would just alter the normalisation B. 

We have performed a least-squares fit to the data. Figure shows the best SM fit which 
corresponds to a mass Mx = 10^^ GeV, with = 132 for 81 degrees of freedom. The low energy 
component (^) and the total spectrum are also shown. In calculating we ignore the highest 
energy point at ~ 3.2 x 10^'^ eV which is based on a single event, hence has a large flux uncertainty. 
Nevertheless this event is important since it requires us to reject values of Mx much below 10^^ GeV 
which cannot generate such an event, although they allow a reasonable fit to the rest of data. In 
Fig. ^ we plot the best SUSY fit taking a common mass of Msusy = 400 GeV for all spartons; the 
favoured decaying particle mass is now Mx = 5 X 10^2 QgY^ ^ith a shghtly lower of 130. 

The assumption of 2-body decay may be rather naive for a superheavy particle like a crypton 
which is expected to decay through very high-order non-renormalisable operators |^^. Many-body 
decay distributes the total energy Mx among several particles and thus fiattens the spectrum. We 
assume that many-body effects are purely kinematical and hence can be encapsulated in the phase 
space of the decay. In Appendix ^ we calculate Pniz), the probability density that one parton 
carries off a fraction z of the total available energy per parton Mx/2. For n > 3 we get 

Pn{z) OC Z{1 - (37) 

To leading order in a<j the particle fiux is then given by 

(z,Ml). (38) 



In particular if D^{x, M^) oc (1 — 1, the differential particle fiux decreases as 

jhalo(^) OC (1 _x)«(*4)+'^-2. 

For many-body X decays the total flux is 

E'J(E) = 4^ + [ '^Pn (f ) D'^-^y^-iz^Ml), (39) 

where n is the number of partons into which X decays. In Fig. ^ we plot SUSY spectra with 
Mx = 10-*^^ GeV, Msusy = 400 GeV and n = 2,8, 16. As n increases the spectrum flattens since 
the average momentum is pushed to lower values of x (see Appendix ^). For n = 2, 8, 16 we get 
= 133, 127, 126 respectively. Clearly many-body decays fit the data better which is encouraging 



as this would be natural in the crypton model |64|. With the forthcoming improvement in event 



statistics a joint fit to the decaying particle mass and decay multiplicity would be warranted. 
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V. DISCUSSION 



We now compare these results to our previous work in which the fragmentation spectra were 
found by directly running HERWIG at the decaying particle mass scale j^^, as well as to other 
work where numerical solutions to the DGLAP evolution were also reported and, finally, 

to a recently proposed independent method for calculating the fragmentation spectrum p^ . 

First it should be emphasised that most studies of 'top-down' models have employed the MLLA 
approximation to the fragmentation functions |^, 45, 45, 91] which, as noted earlier, is valid only 
at X ^ 1 where a very small (and uncertain) amount of energy is released, and is thus incapable 
of being normalised as was done in these papers. Other studies |42| have used a convenient (and 
properly normalised) parameterisation of the fragmentation function based on PETRA data at 
~ 20GeV which, however, does not account for the large scaling violations in evolving up to 
the very high energies of interest here. Moreover in both these approaches it was assumed that the 
ratio of pions to baryons in the cascade remains fixed at its low energy value of 20:1 at all energies. 
(In reality one would expect the FF of say neutrinos to evolve differently from that of baryons due 
to new processes at high energies e.g. decays of heavy flavours.) All these studies focussed on a 
decaying particle mass of ~ 10^^~^^ GeV, close to the GUT scale. 

These difficulties were emphasised in our previous work where HERWIG was used to make 
the first realistic estimate of the fragmentation function for super-massive particle decay, including 
the expected large scaling violations. It was shown that the UHECR data favoured a decaying 
particle mass of 0(10^^) GeV — close to the 'hidden sector' scale of SUSY breaking rather than 
to the GUT scale. However this work suffered itself from the fact that HERWIG overproduces 
baryons by a factor of ~ 3 at high x (see Fig. ||); moreover it did not account for the effects of 
supersymmetry on the development of the cascade. Both these issues were addressed in Ref. ||8^ 
which suggested that the way forward was to evolve measured fragmentation functions at LEP up 
to the decaying particle mass. For the SUSY case, sparton FFs were generated at 10^ GeV using 
the event generator PYTHIA [Q| and the evolution done for a supergravity model with parameters 
Mo = 800 GeV, M1/2 = 200 GeV, ^0 = 0, tan/3 = 10, sgn(^)=+. Apart from highlighting the 
high-x baryon overproduction problem in the previous work for the SM case, this work showed 
that the effects of supersymmetry could be important in that the most likely decaying particle mass 
increased from ~ lO^^ Q^y foj, ^he SM to ~ 10^^ Q^y fo^. gusY Q. Subsequently the DGLAP 
equations were also solved in Ref. ||9^ following a similar procedure but these authors favoured 
the possibility that the decaying particles are uniformly distributed in intergalactic space rather 
than being clustered in the halo. Then the decay spectrum undergoes GZK processing, developing 
a 'bump' below the putative cutoff; by fitting this to the (unnormalised) data from AGASA, Fly's 
Eye and Haverah Park they found an even higher decaying particle mass of 10^^ GeV |^^. (For 
the physically better motivated case of a clustered halo population, their best-fit mass is 10^^ GeV 
in agreement with previous work |4£,86|, but since the FFs for this mass were not presented, we 
cannot compare their results with ours in detail here.) Finally in Ref. |96| a new Monte Carlo 
simulation for jet fragmentation was presented and used to calculate the SUSY FFs for a 10^^ GeV 
mass particle; however the SM FFs were not presented. 
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The authors acknowledge errors in their previous analytic calculation of the SUSY FF in the 
MLLA approximation |91], which was used in several studies of TD models [43,45|. 
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In Fig. we compare the (flavour-averaged) FFs of baryons, photons and neutrinos for a 
10^^ GeV mass decaying particle obtained from HERWIG |4|], with the results from (SM) DGLAP 
evolution in Ref. |8^, as well as in the present study. The baryon overproduction problem at high 
X is evident in the HERWIG results; in fact photons and neutrinos now dominate at all energies. 
Secondly the results obtained by DGLAP evolution agree reasonably well, keeping in mind that 
Ref. 1 86 1 incorporated gluon coherence [76| and also next-to- leading order corrections for both as 
evolution and the splitting functions, which we have not done in the present work. 

In Fig. |ll] we compare the (flavour-averaged) FFs for baryons, photons and neutrinos obtained 
from (SUSY) DGLAP evolution in Ref. |86| and in present work, with the results from the new 
fragmentation code of Ref. |^^. There are significant diflferences between the results. Note that 
we have taken sparton FFs to be zero at threshold rather than generating them with PYTHIA 
above threshold as was done in Ref. |8f:]; this is presumably the major reason for the diflference in 
our results and reflects the present degree of theoretical uncertainty in this approach. The results 
obtained by the new Monte Carlo simulation for jet fragmentation ||9^ diflfer significantly from 
both DGLAP evolution calculations. Since the SM FFs were not presented in this work we cannot 
check whether the new code can reproduce the (consistent) results for SM DGLAP evolution 
discussed above. In view of these discrepancies, further work is clearly necessary to calculate the 
spectrum reliably in the super symmetric case. 

In summary we have developed a new numerical code to solve the DGLAP evolution equations 
for fragmentation functions in heavy particle decay. We have used this to calculate the expected 
spectra of baryons, photons and neutrinos from the decays of hypothetical metastable dark matter 
particles of mass 10^^ GeV clustered in the halo of our Galaxy, both for the case of Standard Model 
evolution and for a simple super symmetric model. The shape of the fragmentation spectrum (of 
either baryons or photons) fits rather well the new component of ultra-high energy cosmic rays 
extending beyond the GZK energy; the fit improves if the effects of supersymmetry and, more 
importantly, of many-body decays are taken into account. 

It is sometimes stated that the decaying particle model is already in conflict with observations 
because the dominant particles in the fragmentation cascades are predicted to be photons rather 
than baryons which are favoured by experiment |19|. However as mentioned earlier the photons may 
be significantly attenuated in their passage to Earth due to interactions on the radio background in 
the halo, particularly if there is a fossil radio 'ghost' around our Galaxy |^^. Whether this is indeed 
possible, subject to the EGRET bound on the low energy 7-rays which will be created by such 
electromagnetic cascades, will be discussed by us elsewhere. Also of interest is the expected flux 
of high energy neutrinos which can be detected in forthcoming experiments such as ANTARES 



and ICECUBE [99|. Most importantly the Pierre Auger Observatory |^2| will significantly 
increase the statistics of UHECRs enabling a better measurement of the spectrum, as well as of the 
composition and (an)isotropy. It will soon be possible to definitively test the exciting possibility 
that cosmic rays have already shown us a glimpse of physics far beyond the Standard Model. 
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APPENDIX A: MANY-BODY DECAY KINEMATICS 



We briefly review how to calculate the phase space Rn for the decay of X into n partons (see 
Ref. 1100(1 for further details); this is needed to calculate the probability density of obtaining a 
decay parton with energy zMxl'2-. 

From a kinematical point of view the decay X — > oi + 02 + 03 + . . . + a„ can be decomposed 
recursively as a 2-body decay X ^ A\ where A\ is a particle with 4-momentum equal to the 
sum of 4- momenta of particles 02, 03 ... a^. Then A\ is decayed into two particles ^1 — 02 + 
where A2 has as 4-momentum the sum of the 4- momenta of 03 . . . a„, and so on up to the decay 
An-2 o-n-i + cLn- This recursivc procedure allows the phase space integral i?„(Mx, mi, . . . m„), 
n > 3, to be written as 

r{Mx-mi)^ dM? riMi-niif (^^2 

Rn= -—^R2{Mx,m,,Mi) ^R^(Mi,m2,M2) 

HA/„_3-m„_2)2 dM2_2 
•••/ IflT R2{Mn-3,mn-2,Mn^2) 'X R2{Mn-2,mn-l,mn), (Al) 

where Mj is the invariant mass of Ai and R2{ma-, mf,, rric) is the phase space integral for the 2-body 
decay a — > 6 + c: 

R2 = ^^X'/\mlmlml)Jdn, (A2) 
\{x, y, z) = + + — 2xy — 2yz — 2xz. (A3) 



We shall study massless partons for which Eq. (Al) reduces to 

f (M| - M?) / |-(M? - M^) 

Mf ^ ^ Jo M| ^ ^ 2^ 

l^(^'-3-^«-2)M^2. (A4) 

Let Pn{z) be the probability density to get parton Oj with an energy fraction z = 2Ei/Mx {Mx/2 
is the maximum energy that a parton can carry, whatever n is). By symmetry this probability is 
independent of the parton considered. Let us choose parton oi. Then z = 1 — M^/M'^ and 

^ (M|- - Mf^ ' 2 ^"-^2 "'^^ 

dM2 



Pniz)dz oc ^(Mi - M/) ^(Mt - M, 



(M^3-M2_2)M2_2- (A5) 



ML2 



The multiple integration can easily be done to obtain 

Pn{z)dz oc (M|- - Mf)(Mf )"-3dMf oc z(l - z)'"-3dz. (A6) 
After normalisation we finally obtain: 

P2{z) = 6{1 - z), (A7) 
p„(z) = (n- l)(n-2)z(l -z)"-^ n > 3. (A8) 

We plot this distribution in Fig. |l^. As the number of final particles n increases, the total energy 
Mx has to be shared between more particles and consequently less energy goes to each one of 
them. The probability density peaks at smaller values of 2; as n grows, accordingly the mean value 
of z decreases as {z)^ = 2/n. 
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FIG. 1. The fragmentation spectrum of protons at the mass peak, as measured by ALEPH, 
DELPHI and OPAL at LEP; the dashed hne is a parametric fit ||8^. The dotted hne is the spectrum 
simulated with HERWIG, illustrating its tendency to overproduce baryons at high x. 
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FIG. 3. Standard Model fragmentation functions for baryons from quarks and gluons, at the 
initial scale Mz (solid lines) and evolved to a decaying particle mass scale of lO^'^ GeV (dashed 
line), 10^^ GeV (dotted line) and 10^^ GeV (dot-dashed line), illustrating scaling violations 
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FIG. 4. Fragmentation functions for baryons from quarks and gluons, at the initial scale Mz 
(solid lines) and evolved to a decaying particle mass scale of 10^^ GeV, for SM evolution (dotted 
lines), and, the more pronounced, SUSY evolution (dashed lines) taking MsusY = 400 GeV. 
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FIG. 5. Dependence of (s)parton fragmentation functions evolved from Mz up to 
Mx = 10^2 Qgv on MsusY = 200 GeV (dashed lines), 400 GeV (solid lines), 1 TeV (dotted lines). 
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FIG. 6. Fragmentation functions for baryons (solid lines), photons (dotted lines) and neutrinos 
(dashed lines) evolved from Mz up to Mx = 10^^ GeV for the SM (top panel) and for SUSY with 
-^SUSY = 400 GeV (bottom panel). 
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FIG. 7. The best SM evolution fit to the cosmic ray data with a decaying particle mass of 



10 GeV. The dotted line indicates the extrapolation of the power-law component from lower 
energies, while the dashed line shows the decay spectrum; the solid line is their sum. 
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FIG. 8. The best SUSY evolution fit to the cosmic ray data with a decaying particle mass of 
5 X 10^^ GeV. The dotted line indicates the extrapolation of the power-law component from lower 
energies, while the dashed line shows the decay spectrum; the solid line is their sum. 
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FIG. 9. Cosmic ray data compared with SUSY evolved spectra for a decaying particle mass of 
10^^ GeV with -MsusY = 400 GeV, and many-body decays to n partons: n = 2 (solid line), n = 8 
(dotted line) and n = 16 (dashed line). 
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FIG. 10. Comparison of fragmentation functions of baryons (top panel), photons (middle panel) 
and neutrinos (bottom panel) obtained from (SM) DGLAP evolution in the present work (solid 
lines) and in Ref. p| (dashed hues), and using HERWIG pi (dotted lines). 
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FIG. 11. Comparison of fragmentation functions of baryons (top panel), photons (middle panel) 
and neutrinos (bottom panel) obtained from (SUSY) DGLAP evolution in the present work (solid 
lines) and in Ref. [p6| (dashed lines), and using a new fragmentation code (dotted lines). 
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